% 参数定义
c = 3.587;
u = 0.25;

% % 实验数据
% T_exp = [55, 55, 100, 25, 25];  % T(P_1), T(P_2), T(P_3), T(P_4), T(P_5)
% r_exp = [2.322, 0, 3, 100, 0];  % r的实验数据
% z_exp = [0, 5.288, 0, 0, 100];  % z的实验数据
% 实验数据
T_exp = [55, 55, 100, 25, 25];  % T(P_1), T(P_2), T(P_3), T(P_4), T(P_5)
r_exp = [2.107, 0, 3, 100, 0];  % r的实验数据
z_exp = [0, 4.2, 0, 0, 100];  % z的实验数据 %

% 定义函数 f(r, z) 和 g(r, z)
f = @(r, z) -abs(z).*log(sqrt(r.^2 + (c + abs(z)).^2) + c + abs(z)) ...
             - abs(z).*log(sqrt(r.^2 + (c - abs(z)).^2) - (c - abs(z))) ...
             + 2*abs(z).*log(sqrt(r.^2 + z.^2) + abs(z)) ...
             + sqrt(r.^2 + (c - z).^2) + sqrt(r.^2 + (c + z).^2) - 2*sqrt(r.^2 + z.^2);

g = @(r, z) log(sqrt(r.^2 + (c + abs(z)).^2) + c + abs(z)) ...
             - log(sqrt(r.^2 + (c - abs(z)).^2) - (c - abs(z)));

% 定义目标函数 (误差函数)
% C_11, C_12, C_2是要优化的参数
objective = @(params) compute_error(params, r_exp, z_exp, T_exp, f, g, u);

% 初始猜测值
initial_params = [1, 1, 1];  % [C_11, C_12, C_2]

% 使用fminsearch进行优化
optimal_params = fminsearch(objective, initial_params);

% 输出训练得到的参数
C_11 = optimal_params(1);
C_12 = optimal_params(2);
C_2 = optimal_params(3);
fprintf('训练得到的参数：C_11 = %.4f, C_12 = %.4f, C_2 = %.4f\n', C_11, C_12, C_2);

% 计算误差函数
function error = compute_error(params, r_exp, z_exp, T_exp, f, g, u)
    C_11 = params(1);
    C_12 = params(2);
    C_2 = params(3);
    
    % 热坏死边界条件的方差项 R_1
    R_1_P1 = C_11*f(r_exp(1), z_exp(1)) + C_12*g(r_exp(1), z_exp(1)) + C_2 + 0 - T_exp(1);  % T0 = 0 for P1
    R_1_P2 = C_11*f(r_exp(2), z_exp(2)) + C_12*g(r_exp(2), z_exp(2)) + C_2 + 0 - T_exp(2);  % T0 = 0 for P2
    
    % 自然条件的方差项 R_2
    R_2_P3 = u*(C_11*f(r_exp(3), z_exp(3)) + C_12*g(r_exp(3), z_exp(3)) + C_2 + 0 - T_exp(3));
    R_2_P4 = u*(C_11*f(r_exp(4), z_exp(4)) + C_12*g(r_exp(4), z_exp(4)) + C_2 + 0 - T_exp(4));
    R_2_P5 = u*(C_11*f(r_exp(5), z_exp(5)) + C_12*g(r_exp(5), z_exp(5)) + C_2 + 0 - T_exp(5));
    
    % 总误差
    error = (R_1_P1^2 + R_1_P2^2) + (R_2_P3^2 + R_2_P4^2 + R_2_P5^2);
end


%预测函数
% % 设定常数和参数
% 
% % 创建 r 和 z 的网格
% r = linspace(0, 20, 200);  % 在 [0, 20] 范围内生成200个点
% z = linspace(0, 20, 200);  % 在 [0, 20] 范围内生成200个点
% [R, Z] = meshgrid(r, z);   % 创建二维网格
% 
% % 计算 evl(r, z) 的值
% Evl = 13.4193 * f(R, Z) -8.9851 * g(R, Z) + 26.6139;
% 
% % 绘制等高图
% figure;
% contour(R, Z, Evl, 20);  % 绘制20条等高线
% colorbar;  % 显示色标
% xlabel('r');  % x轴标签
% ylabel('z');  % y轴标签
% title('等高图：evl(r, z)');  % 图表标题




